//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "statics.hpp"
#include "screen.hpp"
#include "dynamics.hpp"

// Comparing StaticRingIntegrated to StaticRing, trying to find inconsistencies

#define RADIUS 0.33
#define MARGIN_ERROR 2
#define BOX (3*RADIUS)

int main (int argc, char *argv[])
{
  Statics::set_ddsq (0.001);

  Statics stat1;
  ElectroDynamics dyn1(stat1);
  StaticRingIntegrated *pc = new StaticRingIntegrated(1, vect3d(0, 0, 0), RADIUS, 0.001, vect3d(0.0, 0.0, 1.0));
  pc->setCharge (1);
  stat1.addStaticElement (pc);

  Statics stat2;
  ElectroDynamics dyn2(stat2);
  StaticRing *pc2 = new StaticRing(1, vect3d(0, 0, 0), RADIUS, 0.001, vect3d(0.0, 0.0, 1.0));
  pc2->setCharge (1);
  stat2.addStaticElement (pc2);

  prec_t min_perce, max_perce, min_percb, max_percb;
  min_perce = min_percb = 1e9;
  max_perce = max_percb = 0;

  prec_t bstddv, estddv;
  bstddv = estddv = 0;

  int n = 0;
  while (n <= 1000000)
    {
      vect3d p;
      p.x = -BOX + prec_t_drand() * 2 * BOX;
      p.y = -BOX + prec_t_drand() * 2 * BOX;
      p.z = -BOX + prec_t_drand() * 2 * BOX;
      if ( stat1.inContact (p) )
          continue;
      vect3d bf1, bf2;
      vect3d ef1, ef2;
      dyn1.getFields (p, bf1, ef1);
      dyn2.getFields (p, bf2, ef2);
      vect3d bv = bf1 - bf2;
      prec_t perc_b = 100 * bv.length() / bf1.length();
      vect3d ev = ef1 - ef2;
      prec_t perc_e = 100 * ev.length() / ef1.length();

      bstddv += bv.magnitude();
      estddv += ev.magnitude();
      if ( perc_b < min_percb ) min_percb = perc_b;
      if ( perc_b > max_percb ) max_percb = perc_b;
      if ( perc_e < min_perce ) min_perce = perc_e;
      if ( perc_e > max_perce ) max_perce = perc_e;

      n++;
      if ( !(n % 10000 ) )
        {
          bstddv = prec_t_sqrt (bstddv / 10000);
          estddv = prec_t_sqrt (estddv / 10000);
          printf ("Got through %d, min_percb=%g%%, max_percb=%g%%, bstddev=%g, min_perce=%g%%, "
              "max_perce=%g%%, estddev=%g\n", n, prec2double(min_percb), prec2double(max_percb), prec2double(bstddv),
              prec2double(min_perce), prec2double(max_perce), prec2double(estddv));
          min_perce = min_percb = 1e9;
          max_perce = max_percb = 0;
          bstddv = estddv = 0;
        }


//      if ( perc_b < MARGIN_ERROR && perc_e < MARGIN_ERROR ) continue;
//      char b[128];
//      printf ("%s: bfield_err=%g%%, efield_err=%g%%\n", p.sprint(b), (double)perc_b, (double)perc_e);
    }

  return 0;
}

